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Summary 

Improvements have been made to a stream wise upwind algorithm so that it can be used for 
calculating flows with vortices. A calculation is shown of flow over a delta wing at an angle 
of attack. The laminar, thin-layer, Navier-Stokes equations are used for the calculation. The 
results are compared with another upwind method, a central-differencing method, and experi- 
mental data. The present method shows improvements in accuracy and convergence properties. 

Introduction 

Upwind algorithms are important in computational fluid dynamics for calculations of flows 
containing shock waves [1]. Some of them are also able to accurately resolve shear layers 
[2,3]. However, most multidimensional upwind algorithms are first constructed in one dimen- 
sion ami then extended to multidimensions by applying the one-dimensional procedure to each 
coordinate direction. In comparison, the present method uses the local stream direction and 
flow velocity to construct the upwinding. Hence the switching of flux evaluations always takes 
place at sonic values, where the shock waves are located. Therefore this method follows the 
flow physics more closely and in that respect is analogous to the rotated differencing [4] algo- 
rithm developed for the full potential equation. 

The present algorithm is an improvement to a streamwise upwind algorithm, which has been 
applied to steady and unsteady transonic flows over airfoils and wings [5,6,7]. In addition to 
using rotated differencing to implement upwinding in the streamwise direction, the switching 
of fluxes across sonic values is smooth and the entropy condition is automatically imposed in 
a manner similar to Godunov’s method. Contact discontinuities are sharply captured [5] and 
boundary layer profiles are fuller [7] (more accurate) in comparison to a central differencing 
method for a case of separated flow over a wing. 

In the results presented, comparisons are made with the upwind method of Roe [1]. In that 
method, an entropy correction is needed, which results in a convergence difficulty. The present 
method does not exhibit that problem. Other features [7] of the present algorithm are that 
pressure and velocity continuity are enforced in the crossflow direction, and also in the stream- 
wise direction as the velocity approaches zero. These features are adequate for transonic flows 
without the presence of vortices. However for supersonic flows with vortices, two additional 
developments were found to be necessary [8]. First, the manner of using the stream direction 
had to be modified to capture oblique shocks sharply. Second, additional terms were needed to 
stabilize the calculations. The present algorithm is described in detail. 

To demonstrate various capabilities of the present algorithm, the flow was calculated over a 
delta wing with a leading-edge sweep of 75° at a Mach number of 2.8 and at an angle of attack 
of 16°. First, a conical flow approximation was used to limit the calculations to two dimensions, 
so that the influence of grid refinement could be determined. Then full three-dimensional (3D) 
calculations were performed on a medium-density grid. These computed results were compared 
with those of other methods and with experimental data. 



Governing Equations 

The thin-layer Navier-Stokes equations can be written in conservation-law form in a body- 
conforming, curvilinear, coordinate system ( £. rj, C ) as follows: 

d r Q + diE + d v F + d<G = ■^-d< : G v , ( D 

Re 

where Re is the Reynolds number. The vector of conserved quantities Q and the inviscid flux 
vector F are 
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where J is the transformation Jacobian, p is the fluid density, e is total energy per unit volume, 
and H is the total enthalpy. The contravariant velocity component U is defined as 
U = t) x u + 7) y v + riz w. For the £ and £ directions, E and G can be defined similarly. The 
viscous flux vector G v is given in reference [8]. The pressure p is related to the conservative 
flow variables Q through the equation of state for a perfect gas: 

p = (7 — 1) [ e — £(u 2 + v 2 + tt/ 2 )] (2) 

Also, c is the speed of sound, where c 2 = 7 p/p. See reference [8] for the form of equation (1) 
when the conical-flow approximation is imposed. 

Numerical Algorithm 

The upwind algorithm is applied to the inviscid fluxes E, F, and G in equation ( 1 ). (The 
viscous term in equation (1) is discretized by a standard procedure [8], which uses second-order, 
central-differencing.) The upwind algorithm is described by the following formula for the cell 

interface flux F with a surface vector, S = ( rj x , Vz) . 

F(Qi,Q r , S J+ p = x {[ Fi + F r ] + [ F/sign (Ui) + a/ A*F/] cos 2 */ 

1 T 2 J L ( 3 ) 

-[ F-sign (Ur) + s r A*F r ] cos 2 *i- - \A\ AQ sin 2 *} , 

where Qi and Q r are left and right states, respectively, and the metric terms tj x , q„, and rj 2 are 
normalized by |Vtj| as k x = i7 x /|Vt7|, k v = v? y /|Vn|, and k t = tj z /|Vtj|. The contravariant 

velocity U is also normalized by |Vt)| and used as C 7 = k x u + k v v + k z w. For the first-order- 
accurate computations, l = j and r = j + 1 . For higher-order extensions, the MUSCL approach 
[ 3 , 9 ] is used. Sign(C 7 ) equals the sign of U and * is the rotation angle which will be determined 
later. The symbol * indicates local sonic values [ 5 ]. 
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A*F = A*(pq)e a = (p*q* — pg)e* , 
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where q is the velocity magnitude and e s = ( 1 , u, v, w, H) T is the sum of the two acoustic wave 
eigenvectors. Note that equation (4) is based on the rotated difference formula [4,5] for the full 
potential equation. Equation (4) and the switches si and s r , which will be specified, use the 
speed q and the Mach number q/c rather than the velocity component U and the Mach number 
component U/c that many other upwind methods use. With the use of this rotated differencing, 
the switching of terms at transonic shock waves occurs independently of their alignment to grid 
lines. 

The last term in equation (3) is defined as follows: 


\A\AQ =|17|AQ + (c - \U\)[-fe s + pAU e d ] 

A C A (6) 

= -^e, +pcAC/e d + (Ap- =£) |C7| e e + p |C/| e v . 

The variables in equation (6) are averaged between the left and right states, except when 
they follow A . Then, for example, A Q = Q r — Qi. Also, e<j = (0,k x ,k r k z ,U) T , 
which is the difference of the two acoustic wave eigenvectors, e« = ( 1 , u, v, w, q* /2) T and 
e v = (0,e V 2,e W 3,e V 4,e t) 5) :r . Here, e „2 = Au — k x A U, e „3 = Av — k v A U, 
e V 4 = A w — k z A U, and e „5 = ue „2 + t;e „3 + toe V 4 . In Cartesian coordinates, e t is the 
entropy wave eigenvector and e v is a linear combination of the voracity- wave eigenvectors. 

As originally developed [5,6], equation (3) did not use the terms in equation (6). Next the 
terms in equation (6) that use e* and ed were added [7] to enforce pressure continuity in the 
cross-flow direction and in the stream wise direction as the Mach number approaches zero. Fi- 
nally, the terms using e« and e v were added [8] for flows containing vortices. 

Following reference [4], the switches a/ and s r are defined in the manner of Godunov’s 
method as follows: for U > 0 , 


~ 1 — ( 1 — Cm) ( 1 — €r) i 

1 , (7) 

J 1 1+Si g n (^/Tm.r “ 1)1 . 

and M m denotes the Mach number of the averaged state. 

Note that there is current research in improving the sonic point operator [10]. An alternative 
method to those in reference [10] is: for U >0, 

a,= l- €m (2 €i -l) , s r = ( 1 — £ m )( 1 — 2c r ) . (8) 

This smooth switch is identical to equation (7) except at sonic expansion points. At those points, 
one- and two-dimensional calculations using equation (8) have shown increased accuracy over 
equation (7). Equation (8) was derived by modeling a transonic expansion wave for Burger’s 
equation. When the sonic value occurs midway between mesh points, this modeling is exact. 
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Rotated Differencing 

As originally developed, the rotation angle 9 used the cosine of the velocity as cos 6 = U/q 
when SSmlm supeSc. However, i, is importam :to 

to the grid line is beyond the Mach cone. Thus, U/q is replaced by U U/q - U/clfU/c 
become larger than one, cos0 is set to one. This enhances the ability to capture oblique shock 

Wa ^s feature leads to a favorable resolution of bow and crossflow shocks, but it allow j the 
existenceo^CTOSsflow expansion shocks. To avoid expansion shocks, the rotatron angle ts 
determined by a mixture of averaged (m) and pointwise ( l , r) values: 


cos 2 6i r = min[ (1 + 1 1 

Wn W.r 

The following is used for evaluating <f> in this paper because of the smoothness: 

where p, and ps denote upstream and downstream pressures. sin = is deter ‘ 

mined by an arithmetic average of the cosines: sin 0=1- T (cos 6 L + cos 9 r ) . 

Results 

The algorithm given by equation (3) has been tested for flow over a delta wing. Both 2D 
calculation^ 1 usintfthe conic* flow approximation [8], and 3D calculations have been .made, 
compared with B&s method [U] - "“gj 
Computations are carried out in the foUowing manner. The LU-ADI method 1 ” . . 

can be modified for the conical flow fields [13], is used for testing the two upwmd algorithms 
as well as the CD algorithm. Each of the three algorithms is implemented explicitly into th 
LU-ADI method, so that steady flows are determined by each of the thrw algorithms ^mar 
flow is also assumed. For third-order accuracy, the MUSCL scheme with Koren s 
limiter [9] is used. 

Delta- Wing, Conical-Flow Calculations 

This test case considers a vortical flow field over a della wing in order to examine the present 
tatnula's cajmbility for computing *»«--. data 

Figure 1 shows the model geometry, and the typical experimental flow held 
^hnwn schematic allyin figure 2. For the computations, the conical approximanon is used. 

a Refinement study P Tfie coarse, medium, and fine grids all use 

51 poiS normal to the bodf md 27, 51, and 99 pomts < method 
Comoutations were done with the present method. Roe s method, and the CD metnoo. 
Fi^re 3 P shows a comparison of density contour plots of three numerical soluuons on die fine 
SST fThedensity vies are on a portion of the sphere of radius equal to one ) Two shock 
^ ^ i t ho W shock wave on the windward side of the delta wing, 

Wa ;l^, r CT ts :"vfo:L leeward side. The present method and Roe’s 
method giv^ sTrrdUrrcontour plo K for drose sh«k waves, bur .he CD method gives smeared 


(9) 


( 10 ) 
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Dlots (For the CD method, the smoothing coefficient k was set to 0.1 m the fine-grid case 
instead of « = 0 05 in the other cases because at convergence the solution had numerical os- 
ciUations with ° = 0 .05 in the fine grid.) The low density regions at both the primary and 
secondary vortices also indicate that both upwind methods give similar solutions, but the C 

rompSsonof total pressure contour plots on the fine grids is shown in figured 'Hie 
primary vortex appears similar in both upwind solutions, in respect to us location and the con- 
tour level, but not in the CD solution. The primary vortex appears off 
both upwind solutions, but the primary vortex and the boundary layer touch ®* C J} . 

CD solution. The shear-flow region separated from the leading edge also shows diffe :e 
the three solutions. The present formula gives a sharper solution than Roe s in this s ear- ow 

region. Again, the CD formula gives the most dissipative solution. . 

Figure 5 shows comparisons of total pressure profiles normal to the leeward surface o 
wing aty = 0 .216 on the coarse, medium, and fine grids. See reference [8] for comparisons at 
y = 0 1, approximately on the primary vortex, and y = 0 .2, approximately on the secondary 
vortex. Figure 5(c) shows the main features of the complicated profiles. The first peaknear 
z = 0 indicates an edge of the boundary layer under the secondary vortex. The following 
minimum correspond! to the secondary vortex. The next peak near z = 0 .°25 indicates the toud 
pressure recovery between the secondary vortex and the shear layer separated from the leadi g 
edge. This shear layer from the leading edge is observed as the next local minimum. Finally, 
the flow recovers to the free stream. The second peak between the secondary vortex and the 
leading-edge shear layer appears to be higher with respect to both level and location for al 
solutions as the grids are refined. The width between the peak and the region of the recovery 
to the free stream corresponds to the width of the separated shear layer. The present formula 
gives the narrowest shear layer in the three, even on the fine grid. This cnspness indicates that 

the present formula computes the shear flow most accurately. 

Comparisons of pressure coefficient distributions on the leeward wing surface on the coarse 
medium and fine grids aie shown in figure 6. Experimental dam [Ml are also mdicmted m 
figure 6 by upper and lower triangles corresponding to data on the nght- and left-hand side 
of the wing, respectively. Results obtained with the present formula are found to be shg y 
more accurate than those with Roe’s method when compared with experimental data as well as 
with the fine-grid solution. The CD solution has a large discrepancy between the two upwind 
solutions on the coarse grid, but the discrepancy decreases as the grids are refined 

Finally, a comparison of convergence histories of calculations, using the > present _and Roe s 
methods is shown in figure 7. The locally varying time stepping was used [12]. Because ot 
the stiffness of conical source term, A to was set to 0.1 in the first 3000 iterations m e 
impulsive start. Then, At 0 =0 .25 for the next 1000 iterations, and finally At 0 was set to 0.5. 
The maximum CFL number reached about 20. The present formula shows better convergence 
in both the L 2 norm and the J • L. norm (rescaled by the transformation Jacobian) On tire 
other hand. Roe ’s formula reaches a limit cycle. The calculations by both upwmd methods 
started from uniform ffee-stream conditions. However, the calculations by *e CD method 
were started from a converged solution by the present formula. The CD method needed a 
smaller A to : that is, one-fifth of the one used for the upwind computanons. The present formula 
converged the best of the three for this shear-flow computation with resp<*t to the convergence 
rate and the order of magnitude of convergence. The difference in CPU time between 
present formula and Roe’s formula is less than 1% in the present computations. 
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Delta- Wing, Three-Dimensional Calculations 

A medium-density grid in the curvilinear coordinate system ( £, rj, £ ) was used for cal- 
culations for both the present method and Roe’s method. There were 25 x 51 x 41 points in 
the £ (conical), rj (circumferential), and £ (normal) directions, respectively. The convergence 
properties were similar to those in the conical-flow calculations. Both methods produce similar 
results to those obtained in the conical-flow calculations, including the treatment of the bow 
shock wave and the vortices. Figure 8 shows velocity-magnitude plots at an axial location 90% 
of the distance from the nose to the trailing edge of the model. First-order accurate results are 
shown in figure 8(a). The present method (plotted on the right side) shows a slightly larger 
high-speed region than Roe’s method (plotted on the left side). As shown by the third-order 
results in figure 8(b), this higher speed is more accurate. Hence, the present scheme is less 
dissipative than Roe’s. In this 3D test computation, the CPU time per grid point per iteration 
is 36.6, 35.7, and 32.1 fj,sec for the present. Roe’s, and CD computations, respectively, on a 
CRAY X-MP computer. Hence, in 3D, the present method takes about 2.5% more time than 
Roe’s and 14% more time than CD. This confirms that the required arithmetic operations of the 
present formula are comparable to those of Roe and CD. 

Conclusions 

An improved streamwise upwind algorithm has been derived and applied to conical flow 
fields. In comparison with Roe’s method, the present formula (1) captures oblique shock 
waves in the same manner, (2) requires arithmetic operations of the same order, (3) has better 
convergence properties, i.e., no limit cycle, and (4) has an advantage over Roe’s in comput- 
ing shear flows accurately. The results also indicate that the CD method is more dissipative in 
the resolution of shock waves and shear layers than upwind methods. In addition, the present 
method switches differencing at sonic values rather titan at values that are dependent on the 
coordinate system, which is more in accord with the fluid physics. 
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Fig. 1 Model geometry. 


Fig. 2 Experimental flow field 
of vortex with shock 
wave (ref. 14). 
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Fig. 4 Comparison of total pressure con- 
tour plots on the fine grid; 
Moo = 2.8, Re = 3.565 x 10 6 , 
a = 16°. a) the present formula; 
b) Roe’s formula; c) the central- 
difference method. 













Fig. 5 Comparisons of total pressure profiles at y = 0 .216 . a) coarse grid; b) medium grid; 
c) fine grid. 
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Fig. 6 Comparisons of C p distributions on the leeward wing surface, a) coarse grid; 
b) medium grid; c) fine grid. 



Fig. 7 Comparisons of convergence history on the medium grid, a) L 2 norm; b) J ■ Loo norm. 
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Fig. 8 Comparisons of velocity magnitude contour plots, a) first-order; b) third-order. 
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